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Highly acurate numerical simulations are employed to highlight the subtle but important differ¬ 
ences in the mechanical stability of perfect crystalline solids versus amorphous solids. We stress the 
difference between strain values at which the shear modulus vanishes and strain values at which a 
plastic instability insues. The temperature dependence of the yield strain is computed for the two 
types of solids, showing different scaling laws: ~ 7° for crystals versus 7^ ~ 7° _^^y 2/3 

for amorphous solids. 


I. INTRODUCTION 


It is well known that the mechanical stability of bulk 
crystalline solids at finite temperatures is dominated by 
the motion of topological defects like dislocations. In 
perfectly ordered crystalline solids there are no disloca¬ 
tions, and also in amorphous solids the notion of a dis¬ 
location does not exist since there is no long range order 
with respect to which a dislocation can be defined. Both 
crystalline and amorphous solids resist a small external 
stress (or strain) and return to their original shape when 
the stress is removed. On the other hand, when higher 
stresses are applied some brittle solids break while other 
ductile solids exhibit plasticity; they deform and do not 
return to their original shape when the stress is removed. 

Characterizing the mechanical strength of a given solid 
requires an understanding of the values of external stress 
or strain at which the solid becomes mechanically unsta¬ 
ble. We will refer to the values of stress where instabilities 
occur as “critical streses”. For practical purposes one is 
interested in the so-called yield stress (jy which is de¬ 
fined as the highest value of the stress which a solid can 
sustain before undergoing unbounded plastic flow. In a 
generic crystalline solid the yield stress depends on the 
existence of defects, on temperature, on the time of the 
observation, etc. Therefore, in order to define a sharp 
characteristic yield-stress one defines the ideal strength 
- the maximum achievable stress of a defect-free crystal 
at zero temperature. The first attempt to estimate this 
value for an ideal crystal which is elastically unstable was 
made by Frenkel Q, Cf. Eqs. ©-(HSl) below. Recently 
it was shown Q that a crystal can loose stability before 
the critical point predicted by Frenkel, i.e. when one vi¬ 
brational mode reaches zero frequency. In fact, this loss 
of stability occurs before the shear modulus of the crys¬ 
tal vanishes. In this paper we will argue that one major 
consequence of the randomness in amorphous solids is 
that the instability associated with the appearance of a 
soft vibrational mode (zero frequency) is generically after 
the vanishing of the shear modulus. The reasons for this 
important difference will be elucidated and explained in 
Sections ITTTlandUVl 

The critical stresses are calculated at zero temperature 


under quasistatic conditions as is explained in Sect. In] 
In contrast, experiments are usually carried out at finite 
temperatures. Therefore it is important to extend the 
calculation of the critical stresses to finite temperatures. 
In both perfect crystals and amorphous solids the values 
of the critical stresses reduce when the temperature is 
increased, simply because it becomes easier to overcome 
the energy barrier involved in the mechanical instabili¬ 
ties. Nevertheless we will show below, cf. Sect. El that 
the difference between perfect crystals and amorphous 
solids translates to different temperature dependence in 
the reduction of the critical stresses. 

Sect. El presents a summary and conclusions of the 
present paper. 


II. MODELS AND SIMULATION METHODS 
A. Potentials 


In this section we introduce the numerical procedures 
that are common to our analysis of perfect crystals and 
amorphous solids. The different implementations will be 
explained in subsequent sections. 

In all our simulations we employ binary potentials be¬ 
tween pairs of particles. In perfect crystals we have only 
one type of particles, say A, and in the model amorphous 
solids we employ here two types of particles, say A and 
B. The interatomic interactions between particle i (be¬ 
ing A or B) and particle j (being A or B) are defined by 
shifted and smoothed Lennard-Jones potentials 


> \0 if 

( 1 ) 

where 


= 4ey 




( 2 ) 


The parameters are taken from Ref. Q. All the poten¬ 
tials given by Eq. m vanish with two zero derivatives at 
distances i?™* = 2.5aij. The parameters of the smooth¬ 
ing part and details of the interparticle interactions can 
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Subsequently we strain the initial configuration, either 
crystalline or amorphous, with additional external affine 
simple shear. The procedure is as follows: the particle 
positions change under shear strain from the reference 
state {r°} to a new one, denoted {r^}, by an affine trans¬ 
formation that is defined by a matrix J : 

n=J-r°. ( 6 ) 

Here the matrix J in Eq. m is given by J = h( 7 ) • 
h“^( 7 o). It follows from Eq. (|4]) that the matrix J is 
defined by 


FIG. 1: Configuration of the one-component system with per¬ 
fect hexagonal structure. The dotted lines represent the sim¬ 
ulation box which is continued periodically in both directions. 


be found in Ref. [^. It is convenient to introduce re¬ 
duced units, with uaa being the unit of length and caa 
the unit of energy. 


B. The preparation of the initial configuration 

The first step in all simulations is the construction of 
a model solid (crystalline or amorphous) of N particles 
in a two dimensional box of size x Ly with periodic 
boundary conditions. In the case of crystalline solid we 
place the N particles on the vertices of a hexagonal lat¬ 
tice, see for example Fig.[TJ Since the crystal is obviously 
free of defects it is also stress free. Thus the configuration 
is ready for subsequent straining. 

The preparation of the amorphous solid is more in¬ 
volved. Firstly we equilibrate a system with 65% par¬ 
ticles A and 35% particles B at a temperature T = 1 
in Lennard-Jones units. This ratio is chosen to avoid 
crystallization upon cooling. Next we cool the system to 
T = 10“® in steps of AT = until T = and 
then in one step to the final temperature. The obtained 
configuration is not necessarily stress free, with particle 
position denoted by Si from the set Therefore 

we apply simple shear which for a general strain 7 is de¬ 
fined by 

Ti = h{j) ■ Si , (3) 

with the transformation matrix 

h(7)=(jl). (4) 

Note that this transformation is volume preserving. 

The configuration with (almost) zero stress is obtained 
at a strain 70 ; the particle positions at this configuration 
are denoted by {t’°)}(^i, 

= h{-fo) ■ Si. (5) 


where the strain 70 corresponds to the deformation from 
the rectangular simulation box to the reference system. 

In the case of amorphous solid the affine transforma¬ 
tion Eq. ([S]) always destroys the mechanical equilibrium. 
To regain mechanical equilibrium one should allow a non- 
affine atomic-scale relaxation of the particle positions 
{ri} (see, e.g., i). Also for a crystalline solid at fi¬ 
nite temperature one should allow this step of non-affine 
relaxation. At finite temperature this relaxation can be 
performed by Molecular Dynamics or Monte Carlo meth¬ 
ods. In the Monte Carlo protocol one moves the particles 
randomly and the move is accepted with probability 


Ptr = min 


1 , exp 



( 8 ) 


where G is the generalized enthalpy. Under strain control 
the matrix h is fixed and the difference of the general¬ 
ized enthalpy is defined by the difference of the potential 
energy of the system U{h, {a}) 


AG = U{h{j),su--- ,sr^,---SN)- 

U(/r(7),ai,--- ,af",---aiv), 1 < ^ < iV,(9) 


where the displacement of the particle positions is defined 
by 

^new ^ gold l<i<N ( 10 ) 


with the periodic boundary conditions taken into ac¬ 
count. In this equation the a component of the displace¬ 
ment vector of a particle is given by 

= As™,,(2r-1), (11) 

where Asmax is the maximum displacement and is an 
independent random number uniformly distributed be¬ 
tween 0 and 1 . 

It follows from Eq. m and Eq. that in the limit 
T —>■ 0 only the configurations with decreasing energy 
are accepted, i.e., the Monte Carlo process should con¬ 
verge to one configuration with minimal energy. In prac¬ 
tice the direct minimization of the energy of a system 
at zero temperature after every small increase in strain 
(the athermal quasistatic (AQS) strain control protocol 
[a 13) is more effective than the stochastic Monte Carlo 
method. 
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III. HEXAGONAL LATTICE 


A. Thermodynamic instability 

The perfect hexagonal structure is shown in Fig. [T] 
The energy of the system is minimal, U/N = 
—2.5388472, when the distance between neighboring par¬ 
ticles is i?o = 1.12152 (at this point the pressure and the 
internal shear stress are equal to zero) and the dimen¬ 
sionless particle number density is p = 0.918. The de¬ 
pendence of the energy and the shear stress on the shear 
strain 7 under the simple shear defined by Eq. 0 is 
shown in Fig. [51 The elastic shear modulus of the system 
estimated at small strains is equal fi = 24.12. Note that 
the shear modulus vanishes at the maximal and minimal 
points of the stress vs. strain curve in the middle panel 
of Fig. |5J 

The energy is a periodic function of the strain and 
reaches its maximum when the hexagonal lattice is trans¬ 
formed into a square one (which is unstable, see, e.g., Q) 
at the strain 7 = 1/^/{?>). It follows from the stress-strain 
curve (middle panel) that the region between the points 
indicated by square symbols is thermodynamically unsta¬ 
ble. Frenkel proposed an analytical guess for the periodic 
functions shown in Fig. [2] in the form 

rr p(l - cos(\/37r7)) 

^ “ 37r2 

and 

= ■^sm(V37r7) . (13) 

The Frenkel approximation is shown in Fig. |5| by the 
dashed lines. Both the approximation and the numerical 
results indicate that the stress can not exceed some value 
o''‘xy ^ '^Ixy The quantitative details differ. Eq. (USD 
yields the estimation = p/(-\/ 37 r) ~ /j,/ 5, underesti¬ 
mating the result of direct numerical calculation a^y ~ 
/r/4. In fact, Eq. (fT51) and Eq. (fT5]) should be considered 
as first terms in a Fourier expansion Q. The maximum 
value of the stress in the approximation given by Eq. (USD 
corresponds to the inflection point of the strain-energy 
curve at qy = I/(2-\/3) which is associated with theoret¬ 
ical (ideal) strength which is achieved by a homogeneous 
deformation. 


B. Vibrational instability 

1. Pure affine straining 



1 




FIG. 2: The energy (a) and the shear stress (b) under simple 
shear. In blue continuous line we represent the exact, numer¬ 
ically computed data. The dashed red line is the Frenkel 
approximation Eqs. m and (USD. The red triangle and 
the green square represent the vibrational and the thermo¬ 
dynamic instabilities respectively. In panel (c) we show the 
number of negative eigenvalues of the Hessian for the system 
with N = 400 as a function of the strain when nonaffine re¬ 
sponses are suppressed by hand. 


In fact, it is possible to lose stability during purely 
affine straining due to inhomogeneous deformations by 

vibrational modes before becoming thermodynamically eigenvalue of the Hessian matrix. At low temperatures 
unstable. The signifiers of such an instability are the the energy of a system in the solid state can be written 
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FIG. 3: Stress-strain relation for the perfect hexagonal lattice. 
The solid line shows results of AQS simulatins, dotted line 
correspond to the affine transformation (see also Fig. [2|. 


FIG. 4: Lowest eigenvalues of the Hessian for a perfect hexag¬ 
onal lattice with particle number N = 400 and N = 1600 in 
the simulation box. The dashed red lines are an aid to the eye 
to observe the linearity of the dependence of the eigenvalue 
on the strain. 


in the harmonic approximation 


U = Uo + AtTH^At^ , (14) 

with repeated indices summed upon and a, 13 denoting 
the cartesian components. Here Uq is the energy of a 
system in equilibrium and the Hessian H is the matrix 


^ d^U 

dr°‘dr^ 

In a canonical form Eq. (HI reads 


(15) 


For the perfect crystal without defects we expect the 
Hessian to be an analytic function of 7 at least until the 
point of instability. In other words, we can write 


< >= Ap = H( 7 p - 7 ) -h H( 7 p - 7)2 + ... , 

(17) 

where ’®'p is the eigenfunction of the Hessian associated 
with the eigenvalue Ap that vanishes when 7 —^ 7 p. The 
consequences of this analyticity assumption are explored 
below. 


U = Uo + Y.^,Sl (16) 

i 


2 . Relaxational effects 


where Xi are eigenvalues of the Hessian and Si are nor¬ 
mal coordinates. It follows from Eq. HU that in the 
harmonic approximation a solid can be expressed as a 
number of uncoupled oscillators. The structure is stable 
for arbitrary Si if all eigenvalues are positive. The un¬ 
stable deformation begins when the smallest eigenvalue 
approaches zero [lOl - H^ . 

The first eigenvalue Ap crosses zero before the shear 
modulus vanishes, at the value of strain yp denoted with 
the red triangle in Fig. [2] Note that when the strain 
increases this eigenvalue becomes negative, and other 
eigenvalues cross zero and add up to a group of negative 
eigenvalues. The dependence of the number of the nega¬ 
tive eigenvalues on the strain under affine transformation 
is shown in Fig. [5] lower panel. The hexagonal lattice 
loses its stability as a harmonic system much before the 
loss of thermodynamic stability. The reader should note 
that in practice one would never observe this increase 
in the number of negative eigenvalues since the system 
will respond to the instability with non-affine responses 
that are studied next. Here such non-affine effects were 
suppressed by hand. 


The picture obtained with purely affine straining is in¬ 
complete. For more precise and detailed information it 
is necessary to take into account relaxational effects in 
which the system responds to the vanishing of an eigen¬ 
value with non-affine motion. To this aim we apply to 
the same crystalline hexagonal solid an athermal quasi¬ 
static protocol in which after every increase Ay in the 
affine strain we follow up with gradient energy minimiza¬ 
tion to regain mechanical equilibrium [l^ . 

The strain-stress relation obtained in the frame of the 
AQS protocol is shown in Fig. [H One sees that the sys¬ 
tem loses stability before the point of the homogeneous 
instability. It is useful to follow the trajectory of the 
lowest eigenvalue of the Hessian matrix as the strain is 
increased. This is shown in Fig. |4] for two system sizes 
with N = 400 and N = 1600. The point at which the 
eigenvalue vanishes is the same for two system sizes. Near 
this instability point the dependence of A on 7 is well rep¬ 
resented by a linear law. This linearity is a direct conse¬ 
quence of the analyticity assumption (1171) . This will be 
shown to be in marked difference from the amorphous 
solid case. 
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FIG. 6: Instantaneous values of the internal shear stress under 
FIG. 5: Unharmonic model as given by Eq. ([18|). The green strain control for different values of the applied strain, 
triangles denote the extrema of the potential. 


When the harmonic approximation is being lost it is 
necessary to take into account effects of anharmonicity 
in modelling the energy. The simplest model of an an- 
harmonic well is given by 

U{s) = lxph)S^ + (18) 

where Xp{'y) is the lowest eigenvalue of the Hessian and 
K is the constant of the anharmonicity. The dependence 
of the energy given by Eq. on the variable S for 
different Ap( 7 ) is shown in Fig. [5] 

It follows from Eq. (see also Fig. [5]) that the po¬ 
tential barrier is related to the eigenvalue by 


a;7(7) 


2 Ap(7)^ 

3 


(19) 


One should note that Eq. dn is only approximate, 
taking into account only the most unstable mode. In re¬ 
ality, especially in the thermodynamic limit, we expect 
other modes to intervene and dress the predictions dis¬ 
cussed above. This can be seen for example from the 
fact that the first instability shown in Fig. |3] occurs at 
7 Ri 0.15. On the other hand the eigenvalue Ap goes 
to zero at 7 Ri 0.14. Due to the intervention of other 
modes the eigenvalue should become “slightly negative” 
before stability is actually lost. To understand this fur¬ 
ther consider Eq. (HU). Upon the energy minimization 
after the affine step all eigenvalues are effected, some of 
them increase and some decrease. The positive ones add 
to Eq. (ITC)) positively and defer the actual instability. If 
the energy minimization were performed precisely along 
the critical eigenfunction of the Hessian this slight dis¬ 
crepancy would disappear. 


C. Monte Carlo studies at finite temperature 

Monte Carlo simulations are done at finite tempera¬ 
ture, be it as small as it may. This blurs to some ex¬ 
tent the definition of the critical strains associated with 
the instabilities, since temperature fluctuation assist in 
crossing potential barrier. Thus all the critical values 
discussed in this section should be understood as upper 
bounds. It is always possible that longer Monte Carlo 
runs can result in lower value of the critical strains. 

Instantaneous values of the internal shear stress un¬ 
der strain control Monte Carlo simulations are shown in 
Fig. [6l For values of the strain less than some critical 
value the stress fluctuates near a given average value. 
For some critical value of the strain the system dwells for 
some time in a metastable state and then loses stability, 
transforming to a new stable state. We chose the critical 
value of the strain corresponding to the appearance of 
metastable states. 

Results of the Monte Carlo protocol for the mean val¬ 
ues of the energy and shear stress are shown for the crys¬ 
tal in Fig. [T] Under strain control the system undergoes 
a series of transitions associated with a loss of stability. 
Along each elastic branch the system follows the affine 
transformation (with the strain increased by some value 
7 — 7 p), see Fig. [2). Each elastic branch is ending by a 
drop at different values of the strain but with the same 
value of the energy and stress. This values indicate the 
limit of the stability of the hexagonal lattice. With in¬ 
creasing temperature the critical strains decrease. 

At finite temperatures the barrier can be overcome if 
T ^ AU, therefore, the critical value of the eigenvalue is 
given by 

XP{1P) - . ( 20 ) 

The dependence of the lowest eigenvalue of the Hessian 
(for two system sizes) on the strain estimated in the frame 
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FIG. 7: The Monte Carlo results for the energy (upper panel) 
and the shear stress (lower panel) dependence on the strain 
for different temperatures. Circles correspond to simulations 
at T = 10“®, squares to T = 10“^ and triangles to T = 10~^. 



FIG. 8: Temperature dependence of the critical value of the 
strain for the perfect crystal. 

of AQS is shown in Fig. 21 The consequence of the ana- 
lyticity assumption Eq. (1171) is that in the vicinity of the 
point 7 p defined by Xp{'yp) = 0 this dependence can be 
approximated by the linear function Ap( 7 ) = A{^p — 7 ). 
Substitution of this expression to Eq. (EOl) yields 

7 , 7 ° - Ciri/3. ( 21 ) 



FIG. 9: Voronoi diagram for a glass configuration. The color 
code is green for pentagons, white for hexagons and majenta 
for heptagons. Sometime an edge in the Voronoi cell can hard 
to visualize at the scale of this image. 

Results of Monte Carlo indicate the correctness of this 
assessment (see Fig. 0 ). 


IV. MODEL GLASS 

A composition of A and B particles that is stable in 
two-dimensions against crystallization is chosen to be 
65% of particles A and 35% of particles B 0 . The struc¬ 
ture of the configuration of the binary mixture which 
produces our model glass is shown in Fig. [HI 

The typical stress-strain relation of the model glass 
calculated in the frame of the AQS method is shown in 
Fig. [TUI In contrast to the hexagonal lattice (see Fig. [3]) 
instabilities are now appearing at different values of the 
stress. This results from the fact that the hexagonal lat¬ 
tice has only one reference state, in the glass there are 
many reference states and the transition between them is 
caused by a saddle-node bifurcation that is accompanied 
by a sudden drop in stress. 

The fine structure of the stress-strain relation in the 
vicinity of the end of an elastic branch is shown in Fig. 1111 
One can see that there are two special points. One of 
them corresponds to the vanishing of the elastic modulus 
followed by the instability point where the lowest eigen¬ 
value of the Hessian goes to zero. It was shown in @ 
that the lowest eigenvalue of the Hessian tends to zero as 
Ap ^ y/^p — 7 , where 7 p denotes the value of the strain 
at the instability point. When the system is not too 
large and the lowest eigenvalue is well separated from the 
larger eigenvalues of the Hessian matrix it follows from 
this result (which is supported by the simulations) that 
the elastic modulus in the critical region is approximated 
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FIG. 10: AQS stress-strain relation for a glass. The serrated 
line corresponds to AQS simulations with non-affine correc¬ 
tions, the dotted line shows stress-strain relation for a purely 
affine transformation (without non-afhne corrections; the first 
points of instability is indicated by triangles. 

by 

/i Ri /is- , , (22) 

V7P - 7 

where /rs is the Born term. It follows from Eq. (l22ll 
that a theory for the glassy state in the spirit of the 
Frenkel approach would employ for the stress an analytic 
function in the variable x = \/^p — 7 . If applicable, the 
dependence of the stress on strain could be expanded in 
Taylor expansion around the point 7 p 0 

oxyin) = <^p + Yl 

i=l 

where C 2 = /ip. In fact this expansion may not exist 
and higher order term may diverge in the thermody¬ 
namic limit due to the accumulation of small eigenvalues 
of the Hessian (prevalence of many low lying barriers), 
as demonstrated in Ref. [ 2 ^ 


A. The difference between crystal and glass 

Both for the hexagonal lattice and the glass there is a 
point of instability defined by a vanishing shear elastic 
modulus (point A). Another instability point (point B), 
related to vanishing the lowest eigenvalue of the Hessian 
appears before point A in the stress-strain dependence 
of the hexagonal lattice but after point A in the case of 
glass. This difference has the following consequence: in 
the case of the hexagonal lattice when the strain is lower 
than point A the system is thermodynamically stable, 
and there will be no important difference between stress- 
controlled and strain-controlled protocols. In both cases 
the stress can be equilibrated in the system such that 



T 


FIG. 11: The shear stress (upper panel) and lowest eigen¬ 
value of the Hessian (bottom panel) dependence on the ap¬ 
plied strain for a glass configuration. Note that in this case the 
point A (denoted by the triangle) where the shear modulus 
vanishes precedes point B where the Hessian lowest eigenvalue 
\p goes to zero. 
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FIG. 12: Dependence of the critical strain value on the tem¬ 
perature for a glass. 


in stress-controlled protocols the internal and the exter¬ 
nal stress are equal. Accordingly one can expect a sim¬ 
ilar temperature dependencies for 7 ^ (T) under stress or 
strain control. 

In contrast, in a glass under stress-control protocols 
the vanishing of the shear modulus is defined by point A 
with the lowest eigenvalue of the Hessian being still finite. 
Therefore, imagine that we apply to the glass a stress- 
controlled protocol with the external stress being smaller 
than the critical stress at point A. At this situation the 
systems is still experiencing a barrier that needs to be 
overcome since \p Q. AtT = 0 therefore we will not 
experience an instability. 

The temperature dependence of the strain critical 
value obtained in the frame of the Monte Carlo proto¬ 
col is shown in Fig. 1121 The temperature dependence of 
the yield strain is in agreement with ~ t 2/3 behavior 
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FIG. 13: The dependence of the lowest eigenvalue of the Hes¬ 
sian on the applied strain for a glass configuration under affine 
transformation. 

V. CONCLUSION 

We have presented highly accurate numerical simula¬ 
tions to underline some fundamental difference between 


the instabilities of glassy materials and perfect crystals, 
even when the atomistic interaction are the same. The re¬ 
sults indicate the importance of examining small systems 
where the precise profiles of the stress vs. strain curves 
can be visualized. Increasing the system size results in 
reducing the strain or stress differences between points of 
instability, and eventually obliterating the details of the 
precise form of the stress vs strain characteristics. 

Fundamentally, the difference is in the analytical de¬ 
pendence of the eigenvalues of the Hessian matrix on the 
strain (or the stress). We note for example Fig. [TUI where 
we highlight the distinction between straining the system 
allowing non affine response and not allowing it. In the 
first case the eigenvalue has a square-root singularity as 
a function of the strain, as discussed in Sect. IIVI In the 
second case, cf. the dotted linear in Fig. [TUI the lowest 
eigenvalue of the Hessian matrix vanishes in an analytic 
fashion, liner in the strain, much in the same way as in 
the crystalline case, cf. Fig [TUI The avoidance (by hand) 
of the saddle-node instability of the non affine response 
results in a fundamental change in the analytics of the 
dependence of the stress on the strain. 
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